Josephson current in graphene: the role of unconventional pairing symmetries 
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We investigate the Josephson current in a graphene superconductor/normal/superconductor junction, where su- 
perconductivity is induced by means of the proximity effect from external contacts. We take into account the 
possibility of anisotropic pairing by also including singlet nearest-neighbor interactions, and investigate how 
the transport properties are affected by the symmetry of the superconducting order parameter. This corresponds 
to an extension of the usual on-site interaction assumption, which yields an isotropic s-wave order parameter 
near the Dirac points. Here, we employ a full numerical solution as well as an analytical treatment, and show 
how the proximity effect may induce exotic types of superconducting states near the Dirac points, e.g. p x - 
and p H -wave pairing or a combination of s-wave and p + lp-wave pairing. We find that the Josephson current 
exhibits a weakly-damped, oscillatory dependence on the length of the junction when the graphene sheet is 
strongly doped. The analytical and numerical treatments are found to agree well with each other in the s-wave 
case when calculating the critical current and current-phase relationship. For the scenarios with anisotropic 
superconducting pairing, there is a deviation between the two treatments, especially for the effective p^-wave 
order parameter near the Dirac cones which features zero-energy states at the interfaces. This indicates that 
a numerical, self-consistent approach becomes necessary when treating anisotropic superconducting pairing in 
graphene. 



PACS numbers: 74.20.Rp, 74.50.+r, 74.20.-z 

I. INTRODUCTION 

The unusual electronic properties of the charge-carriers in 
graphene^ have triggered a massive interest in this material 
over the last few years. Graphene is a monolayer of graphite, 
and thus has a two-dimensional honeycomb lattice structure 
consisting of two triangular sublattices. The two most inter- 
esting features of the dispersion relation for the quasiparti- 
cles moving in a graphene sheet is that (i) the quasiparticles 
at Fermi level are nodal, meaning there is no Fermi surface at 
zero doping, and that (ii) the band structure is conical, thus 
giving rise to an effective mass of zero for the quasiparti- 
cles. These facts have paramount implications for a number 
of physical properties of a graphene sheetP^l 

Quite recently, proximity-induced superconductiv- 
ity in graphene was achieved experimentally by means 
of d epos iting superconducting contacts on a graphene 
snee (|4|5|6j rp n i s nas i e( j t multiple investigations with 
respect to the transport propertie s of superconducting 
grapheneP M ' 11 ' 1 ^ 13 ' 14 ' 15 ' 16 ' 17 ™^ Some of the key 
findings in these investigations include the possibility of 
specular Andreev reflection^ oscillations in th e c onduc- 
tance of a superconductor/normal (SN) junction, 9 1 1 and a 
finite Josephson current even at the Dirac point in an SNS 
junction. 8 Very recently, some authors have also expl ored 
the interplay between proximity-ind uced ferro magnetisrrl^^ 
and superconductivity in graphene j 23 l 24 l 25 l 26 l So far in the 
literature, conventional s-wave superconducting pairing has 
been the primary focus, whereas only little attention has 
been paid to how unconventional pairing in superconducting 



graphene structures influence the transport properties ! 15 ! 19 ! 

In the majority of studies considering transport properties 
of superconducting graphene hybrid structures only an ana- 
lytical scattering matrix approach has been employed. The 
advantages of such a treatment as compared to a purely nu- 
merical one is that it often offers more physical insight into 
the problem. On the other hand, a numerical self-consistent 
treatment is more authoritative and will in general provide 
more accurate results both qualitatively and quantitatively.^ 
Ideally, it would therefore be desirable to compare an analyti- 
cal approach with a numerical treatment in order to see which 
conclusions obtained in the former case still hold in the latter 
case. 

Motivated by this, we present in this paper both an ana- 
lytical and numerical study of the Josephson current in an 
SNS graphene junction (see Fig. [TJ using both conventional 
and unconventional superconducting contacts. We here allow 
for the unconventional pairing by means of including nearest- 
neighbor interactions in a tight-binding model. This interac- 
tion gives rise to exotic types of superconducting states near 
the Dirac points, which also significantly alters the behavior 
of the supercurrent in the system. Our main result is that, 
whereas the analytical and numerical treatments are found to 
agree well with each other in the isotropic s-wave case, there 
is a deviation between the two treatments in the anisotropic 
case, especially for an effective p^-wave order parameter 
near the Dirac cones featuring zero-energy states. This find- 
ing suggests that a numerical, self-consistent approach is re- 
quired when studying anisotropic superconducting pairing in 
graphene. 

The article is organized as follows. In Sec. |ll] we lay the 
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FIG. 1: (Color online) The experimental setup proposed in this 
paper. Two superconducting electrodes in close proximity to a 
graphene sheet induces superconductivity. The Josephson current 
flows through the graphene sheet indicated by the red arrow. Top and 
bottom gates contacted to the graphene sheet permit local control 
over the chemical potential. 




FIG. 2: (Color online) The graphene honeycomb lattice with the two 
different atomic sites A and B, the three nearest neighbor directions 
{ai, a,2, £13}, and the zigzag and armchair interfaces marked. 



theoretical foundation in terms of notation and formalism. 
Sec. 



In 



III we present and discuss our main results, focusing 



on the analytical treatment in Sec. |III A|and the numerical ap- 



proach in III B Finally, we summarize our findings in Sec. IV 



H. THEORY 

A. Analytical treatment 

Our starting point is the tight-binding Hamiltonian on the 
graphene lattice, including a superconducting pairing order 
parameter which is induced by the proximity effect from ex- 
ternal superconducting contacts. We include the possibility 
for both on-site superconducting pairing and nearest-neighbor 
pairing by means of nearest neighbor spin-singlet bond (SB) 
correlations. The full Hamiltonian takes the form: 

ff=-tj;(4^,,iEc) 

ij 

+£[a £/ (4 t 4i+44)+ h - c -] 



constant terms. 



(1) 



Here, Ai a and Bi„ are the second quantized fermion opera- 
tors on the sublattices A and B, while a,j denotes the three 
nearest-neighbor vectors, see Fig. [2] They read a\ — o(l, 0), 
0,2 = a(— l,v3)/2, 03 = a(— 1, — v3)/2, where a is the 
inter-atomic distance of the carbon atoms. For later use, we 
also define the reciprocal vectors to the Dirac points K± 
as K± = (0, ± 3 ^ a ). The superconducting pairing is ac- 
counted for by Ajj, which corresponds to the on-site inter- 
action, and Aj a ., which corresponds to the pairing interac- 
tion along the nearest-neighbor vectors aj, Aj a , may in gen- 



eral be different for different a/s, i.e. the three different near- 
est neighbor bonds. In fact, a self-consistent solution admits 
three possible solutions 27 classified as follows upon defining 

Aj = (Aj ai ,Aj a2 ,Aj a3 ): 



Extended s-wave: Aj = A t (l, 1, 1), 
d x 2_ y2 -wave: Aj = A t (2, -1,-1), 
d^-wave: Aj = A t (0, 1, -1). 



(2) 



The classification of the different symmetries stems from 
which irreducible representation of the crystal point group D 6 
they belong to when considered in the whole band structure 
Brillouin zone, i.e. the full reciprocal unit cell in the basis 
where the kinetic energy is diagonal. The extended s-wave 
gap is proportional to the band dispersion, i.e. oc e q given 
below in Eq. ( 30 1. It has the full symmetry of the lattice, 
thus belonging to the Ai irreducible representation although, 
in contrast to the on-site s-wave Ay, it varies in magnitude 
over the Fermi surface. The different d-wave solutions have 
four-fold symmetries and belong to the two-dimensional E2 
irreducible representation and, therefore, technically any lin- 
ear combination of these two solutions is a valid solution from 
a symmetry standpoint. For an effective potential giving rise 
to an intrinsic SB pairing in the translational invariant bulk 
the two d-wave states are degenerate at T c but the complex 
combination d x 2_ y 2 + \d xy is favored just below T C W^ Inter- 
estingly, this state breaks thus time-reversal symmetry (TRS). 
Here, however, the SB pairing is induced into the graphene 
from external contacts and we choose to limit the symmetries 
studied to the ones given in Eq. Q as those would be the 
ones naturally induced from correspondingly aligned <i-wave 
superconducting contacts, such as high-T c cuprate supercon- 
ductors. 

Introducing the Fourier transform of the fermion operators 
according to 

Aqa — ^ ' Ai a e q ' , Ai a — S ' A qcr e q 1 , (3) 



and similarly for A — > B, we may write down the Hamilto- 
nian in momentum space (now discarding irrelevant constant 
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terms): 
H 
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V 

1 

( A lv B iv A -ii' B -"i)> 



Af„ = 



-fi 



At(-g) 



.t 



A(g) \ 
A(-g) A y 

-e(-g) A* / 



(4) 



Again, Ay is the superconducting order parameter resulting 
from on-site pairing interaction, while we have defined 



e(q) = A(q) = ^ A Ja .e^ 



(5) 



We are interested in the behaviour near the Dirac points K±, 
and hence wish to evaluate M„ at q = -RT-t + fc where |fe| <C 
|-K±|, i.e. in the low-energy limit. In this case, we find by 
means of a straight-forward Taylor expansion, that 



e(K± + k) = v F (-ik x ± k y ), 



(6) 



with the definition v F = 3ta/2, while for the superconducting 
order parameters we have 



Extended s-wave: A(K± + k) 



3A t a 



i}k x T ky), 



d x 2_ y 2-wave: A(K ± + k) = 3A t [l + -{\k x ± fc y )], 

da-^-wave: A(.K"± + fe) = \/3A t [±i + — (— \k y ± fc x )]. 

(7) 



From Eq. the classification of the different symmetries for 
nearest-neighbor pairing is far from obvious. This is because 
in Eq. |7]i they are expressed in reciprocal space where the 
kinetic energy is not diagonal. By diagonalizing the kinetic 
energy through an unitary transformation on A and B, the 
correct symmetries will appear. These low-energy expansions 
are given in Sec. |III A 2| Diagonalization of the Hamiltonian 
in Eq. Q yields the Dirac Bogoliubov-de Gennes (DBdG) 
equation which describes the quasiparticle excitations in the 
system. We find that the DBdG equation close to the Dirac 
points may be written as 



e(K ± + k) 
fe) -fj, 
v A^K T -k) 
\A\K ± + k) A 



( 

At 



4, 



A(K T - 
M 

-e(K T - 



fe) 
■fe) 



A(fsT±4 
A* 



fe) \ 



tp = Eip. 



(8) 



/ 



Due to the valley degeneracy, it suffices to consider only the Dirac point K+. For concreteness, we include on-site and extended 
s-wave symmetry superconducting pairing below. It should be noted that the superconducting pairing potential couples electron 
and hole excitations between the two valleys K + and K_ in Eq. ([8]), as required by time-reversal symmetry (within a single 
valley, time-reversal symmetry is broken). In this case, we find that 



/ —fi v F (-ik x + k y ) 

V F (lk X + ky) -/i 

A\j 3A\a(ik x - k y )/2 

V-3A|a(ifc x + ky) /2 At 



A v 3A 
-3A t a(ik x + k y )/2 A v 

fi -v F (-ik x + k y ) 
-v F {ik x + k y ) fi ) 



tp = Eip. 



(9) 



Note that normal-state (non-superconducting) contribution to 
the above is slightly different from the usual Dirac equa- 
tion — fiao + v F p ■ cr. In fact, the upper-left 2x2 matrix 
of Eq. (|9|l may be written as — (xctq + v F (k x a y + k y a x ). 
The reason for this discrepancy is that the exact form of the 
Hamiltonian depends on the choice of nearest-neighbor vec- 
tors a, , j 6 {1,2,3}, that were introduced previously. The 
physics must clearly remain completely unchanged regardless 
of the choice of aj. However, to facilitate comparison with 
previous work in the literature we revert to a choice that yields 
a Dirac-like Hamiltonian for the normal-state. This is sim- 
ply accomplished by switching the coordinate system chosen 
originally, 1.6. i y k y . Performing this substitution in Eq. 



(|9]l yields the desired form of the normal-state Hamiltonian. 

The strategy for calculating the Josephson current in the 
junction is to first obtain the energy spectrum for the An- 
dreev bound states in the normal region of graphene. This 
is done by matching the wavefunctions at the two SN inter- 
faces, and then solving for the allowed energy states. Explic- 
itly, the boundary conditions dictate that ^lI^^o = ^n\x=o 
and ^r\ x —l = ^n\x=l, where L is the length of the N re- 
gion, i.e. the junction length and 

v N = hip% + t 2 ipt + t 3 ip\ + t 4 ip b _, 

We allow for the chemical potential to be different in the S 
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and N regions. Finally, note that the subscript ± on the wave- 
functions in the normal region indicates the direction of their 
group velocity, which in general is different from the direction 
of momentum. Consequently, although the Andreev-refiected 
hole wavefunction carries a subscript "— " above, one should 
keep in mind that for normal Andreev reflection, the direction 
of momentum is opposite to the group velocity for the hole. 

The Josephson current is computed via the usual energy- 
current relation summed over projections of all paths perpen- 
dicular to the tunneling barrief^ 



iAW) 



4c 



E 



tt/2 



-tt/2 



d7COS7 dSi(A(j)) 



(11) 



where £j(A<^) are the Andreev bound states carrying the cur- 
rent in the N region, and A(j> = 4>r — 4>l is the macroscopic 
phase difference between the superconductors. The integra- 
tion over angles 7 takes into account all possible trajectories 
and f(x) is the Fermi-Dirac distribution function. We define 
the critical supercurrent as I c = \max{Ij(A<p)}\, and note 
that the factor of 4 in front of the summation in Eq 



11 



is 



due to the spin- valley degeneracy. The formula for the Joseph- 
son current disregards the contribution from supergap states, 
which is allowed as long as L/£ <C 1, i.e. a short junction. 



B. Self-consistent numerical treatment 

A self-consistent numerical treatment allows for spatially 
varying order parameters, Ajj and A,;, and will thus directly 
capture the proximity effect inside the junction through the 
depletion of pair amplitude in the superconductor near the in- 
terface and the induction of pair amplitude into the normal 
region. From this it is also possible to explicitly calculate 
the full Josephson current without any restrictions to small 
junctions as is the limitation for the analytical treatment with 
Andreev bound states. Naturally, the local density of states 
(LDOS) will also be readily available. Here, we will use 
the tight-binding Bogoliubov-de Gennes (TB BdG) formal- 
ism, which allows for a self-consistent solution of the order 
parameters as a function of position. The procedure has been 
outlined in detail in Refs. [12] and [29 1 and we will here only 
outline the essentials in order to connect to the analytical treat- 
ment and interpret the results. 

As in the case of the analytical treatment above of a 
graphene SNS Josephson junction we only want to model the 
actual graphene sheet. Therefore, we have to capture the effect 
of the superconducting contacts deposited on top of the sheet 
by some effective parameters. In the analytical treatment this 
was simply done by assuming that the contacts induce con- 
stant order parameters, or gaps, Ajj and A,/, in the S regions. 
For spatially varying order parameters we need to go beyond 
that approximation and we will instead model the effect of the 
the superconducting contacts by using effective pairing poten- 
tials which are only nonzero in the S regions of the graphene 
sheet. For s-wave contacts the pairing correlations induced 
into the graphene are modeled by a simple attractive Hubbard 
J7-term. The nearest neighbor SB pairing can in the same 



way be produced by an effective SB potential J where the ef- 
fective coupling is given by a JS\ ■ Sj term between nearest 
neighbors.^ Thus, the starting point for a TB BdG treatment 
is the following effective HamiltoniarJ23 



ff eff = -f^(iB i+0 ,„+li.c) 

ijcr 

i 



where 



F; 



t _ 



(13) 



is the nearest neighbor SB creation operator. With a simple 



Hartree-Fock-Bogoliubov mean-field approximation Eq. ( 12 1 
can be transformed into Eq. ([TJ, but now with the position 
dependent order parameters 



A,(i) = -(7(i)Mt)l<W 



(i) = - J(i)<^B i+ajT - AifB i+ajl ) 



(14) 

i+a,\l- (15) 



A standard TB BdG formulation of this mean-field Hamil- 
tonian, Eq. ([lj, will result in a 4A?~ x 4A*" eigenvalue prob- 
lem, where N is the number of unit cells in the whole junc- 
tion, and position-dependent, BCS self-consistency equations 
for Ajj(i) and Aj(i). By starting with an initial guess for 
the order parameters Ajj(i) and A,/(i), then solving the 
eigenvalue problem for these values, and finally using the 
BCS self-consistency equations, we can compute new val- 
ues for the order parameters and continue the process until 
self-consistency is reached. For a specific symmetry of the 
SB pairing contacts, we simply restrict A, 7 to that particular 
symmetry. Since the order parameters are by definition zero 
in the N region where U and J are zero, the relevant param- 
eters to describe the proximity effect inside the junction are 
instead the pairing amplitudes Fjj(i) = —Ajj(i)/U(i) and 
Fj aj (i) = — Aj aj (i)/[v2J(?)]. The other significant quan- 
tity to study in a SNS junction is the Josephson current. We 
calculate this by fixing the phase of the order parameters in 
the very end of the contacts and then solving self-consistently 
in the rest of the sample. The Josephson current can be calcu- 
lated relatively straight-forwardly using the continuity equa- 
tion and the Heisenberg equation for the electron density, see 
Ref. Ifl2l for further details. It should be noted that in this 
approach the phase of the order parameter is allowed to vary 
even in the S regions, except at the very end of the contacts. 
This ensures true bulk-like conditions and gives a consistent 
Josephson current throughout the structure. However, it has 
the side-effect that the phase difference Acf> over the junction 
itself will always be less than it, since that is the largest phase 
difference we can apply across the whole structure, but part 
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of this drop will necessarily take place in S if the current is 
non-zero. While this appears as a numerical artifact in this 
context, it is in fact closely related to the physical 2ir phase- 
slip process in Josephson junctions (see e.g. Ref. I30l0 . 



1. Simulation details 

Since the Aj superconducting state corresponds to pairing 
along the nearest neighbor bonds on the bipartite honeycomb 
lattice there will be a directional dependence for this state. 
In particular, different interfaces will behave differently. The 
two most common interfaces for the honeycomb lattice are the 
zigzag and the armchair interfaces, see Fig. [2] but any chiral 
interface is experimentally possible. Now, since the particular 
symmetries for the Aj state, Eq. (|2b, are either fully sym- 
metric or have a four-fold symmetry, the number of interface 
vs. symmetry combinations can be reduced. Obviously, for 
the extended s-wave the interface orientation does not matter. 
For the d-wave solutions, the solution for d x 2_ y 2 -wave on a 
zigzag interface should, at least to a good approximation, be 
the same as the d^-wave on an armchair interface. But, as 
will be shown, the d x i_ y i- and cL^-waves on the same inter- 
face behave significantly different. It is therefore of interest 
to study both symmetry solutions but it suffices to look at one 
interface. The analytic results in this article are calculated 
for an interface along the y-direction, i.e. the zigzag interface, 
though the detailed shape of the interface is obviously irrel- 
evant in a continuum model. For the numerical treatment, it 
turns out, however, that the d^-wave is strongly suppressed 
at external edges, thus demanding bigger contact regions, so 
instead of studying the d x 2_ y 2- and d xy -waves on the zigzag 
interface, we instead study the 6^2 -wave state on both the 
zigzag and armchair interface, with the latter equivalent to the 
d^y-wave solution on the zigzag interface. In order to reduce 
the computational cost we consider only clean, smooth inter- 
faces and Fourier transform in the direction parallel to the in- 
terface, thus significantly reducing the number of unit cells 
N. From an experimental point of view this entails studying 
junctions with an infinite width, though the periodicity of the 
solution is limited to one unit cell. 

The physical input parameters in the TB BdG treatment 
are the on-site pairing potential U for conventional s-wave 
contacts, the SB pairing potential J for unconventional con- 
tacts, the effective potential ji in S and N, the length L of 
N, and temperature. For conventional s-wave superconduct- 
ing contacts we choose the following setup: U(S) = 3.4 eV 
= 1.36*, n s = 1.5 eV = 0.6t. This leads to A v = 0.1 eV in 
the bulk which corresponds to a superconducting coherence 
length £ = Hvp/Au ss 50 A which is 25 unit cells in the 
zigzag direction and 40 unit cells in the armchair direction. 
These values satisfy Xp(S) <§; £, allow us to numerically in- 
vestigate both the L < £ and L > £ cases, and coincide with 
previous work.El They are, however, quite large values for a 
realistic situation but a smaller superconducting gap leads to 
a slower convergence rate and also a need for a larger system 
making calculations less feasible. We have checked our key 
results for smaller U and found no significant difference. 



In order to be able to compare the results between the con- 
ventional and unconventional contacts we need to use the 
same strength superconductor. This is not completely triv- 
ial since the unconventional contacts have an energy gap that 
varies on the Fermi surface. For the proximity effect in a 
SNS junction the relevant scale is the coherence length £ 
since this parameter will determine the superconducting decay 
length. We therefore choose the effective pairing potentials 
such that £ is unchanged between the different symmetries. 
For the unconventional superconducting contacts we calcu- 
late £ = Kvp/Ap, with Ap = (A|)pg being the average 
of the energy gap Afc over the Fermi surface. 31 Since the dop- 
ing level is very high in the contacts, we can approximate the 
band structure to only consist of one band and the energy gap 
parameter for the unconventional contacts then simplifies tcpZl 

A fc = ^2 A Jo,j cos ( k ■ a 3 - <Pk), (16) 

3 

where <pk — arg(^ ■ e lk aj ). Using this we find A^(d-wave) 

= ^ A j where A j is the norm of Aj. Then finally, by com- 
paring the coherence lengths, we can set J(S) = 2.45 eV for 
a straightforward comparison with the conventional s-wave 
contacts. Please note that in the analytical formalism we need 
to choose A t different for all three symmetries in order to have 
the same £. For example, if Ajj = A then A t (ext. s) = |A , 

A t {d x 2_ y 2) = ^A , and A t (d xy ) = ^§A , respectively. 

We have for simplicity assumed that the doping profile 
changes abruptly from /is to //n at the interface. This is 
the same approximation as in the analytical solution and will 
therefore provide the most accurate comparison. We have 
studied junctions ranging from the undoped regime, /Ltjv = 
eV, to moderately doped, = 0.7 eV, to the case of no 
Fermi wavevector mismatch (FVM), i.e. /i N = 1.5 eV. 

In terms of L, we have studied zigzag junctions with L = 2- 
60 unit cells (1 unit cell = 2.13 A) and armchair junctions with 
the corresponding number of cells (1 unit cell = 1.23 A). The 
temperature was chosen to be T = 10 K throughout the work, 
which in comparison to T c in the S regions is effectively zero 
temperature. The accuracy of the solution is determined by the 
choice of termination criterion for the self-consistency step, 
the number of fc-points used in the Fourier transform, the size 
of S and the region in S where the phase of the order paramter 
is kept constant, all of which have been tested thoroughly. 



III. RESULTS AND DISCUSSION 

A. Analytical results 

In order to investigate how the unconventional pairing near 
the Dirac points affects the Josephson current, we begin by 
obtaining an analytical solution for the most possible gen- 
eral case. In a realistic situation, there could be both on- 
site and nearest-neighbor interactions, thus giving rise to both 
isotropic s-wave pairing and one of the order parameters given 
in Eq. {7}. 
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1. Extended s-wave 

We first consider the extended s-wave case, which gives 
rise to an effective p + rp-wave order parameter near the Dirac 
points. The BdG-equation then reads 



—(j, pe 
pe 10 -n 
Ate" 1 ' 



Ac/ 
A T e ie 

A* 

-pe lB 



A T e- l9 \ 
Ac/ 



-pe 
/' 



(17) 



/ 



We here use units h = vp = 1 and have defined p x ± ip y = 
pe ±l6 and A^ = —3pA t a/2. The above matrix is Hermitian 
as required and may be diagonalized to yield the eigenvalues 



£ a p = aV(/3p-M) 2 + |At/ + /3A T | 2 , a,p = ±. (18) 

Here, j3 = ±1 denotes the conduction and valence band while 
a = ±1 distinguishes between electron-like and hole-like ex- 
citations. The eigenfunctions may then be constructed in the 
superconducting regions. We here only consider positive ex- 
citation energies e > 0, and impose the mean-field restric- 
tion that the chemical potential /is in the superconducting re- 
gions must be much larger than the superconducting gap, i.e. 
/is ^ A. In this case, only the conduction band j3 = 1 
partakes in the low-energy scattering processes. We define 
[v(p)] 2 = 1 - [u(pj\ 2 with 



«(p) = VK 1 



v/e 2 -|Au + A T (p)| : 



(19) 



From Eq. ( 18 I, one can then solve for the wavevector as fol- 
lows (a = 1 since e > 0): 



fJ-s 



A V A T ± 



<?e,h 



A T ) - ((j,At + A[ 



Af, 



(20) 



where A^ = At/p — —3A t a/2 and the scattering angles 
are obtained by means of translational invariance in the y- 
direction: 



qi sin ( 



p e sin ( 



e,h 



(21) 



We also define A/i = /j,g — hn and 



e + fx N , p h 



Pat , p h sin 9 A = p e sin 0, (22) 



related to the quasiparticle momenta in the normal graphene 
region. Here, Afi is the difference in the local chemical poten- 
tial between the S and N region, which may be experimentally 
controlled by means of a gate voltage on top of the normal 
graphene segment. 

Finally, we are able to write down the wavefunctions in 
the three regions shown in Fig. [T] We remind the reader 
that the effective pairing near the Dirac points is a mixture of 
isotropic s-wave and p + lp-wave. However, note that while 
the extended s-wave in this picture corresponds to the even 



and complex p + lp-wave symmetry, it breaks neither time- 
reversal symmetry nor is it a spin-triplet. The deceptive ap- 
pearance is due to the fact that the honeycomb lattice has two 
distinct Fermi surfaces at all doping levels considered here 
and when considered together, the spin-singlet character and 
TRS invariance is preserved as it should be for an extended 
s-wave. Note that in what follows, we will assume that the s- 
and p + lp-wave superconducting order parameters are phase- 
locked, i.e. they are characterized by the same broken U(l) 
gauge phase. For x < 0, we have 



S,L 



+ t h L 



( u(q e ) \ 
u{q e )e^- e ^ 

/ v(q h ) \ 
v(q h )e lSh 
u{q h )e~^ 
\«(?fc)e , ^-* 1 V 



-iq e cos u e x 



lQh COS UfrX 



(23) 



^ N = a 




gip e COS Bx , £ 




ip e COS vx 



(24) 



+ t h R 
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u(q e )e ie " 
v(q e )e-^ 
Vw(g e )e I ^-^«) 

/ <q h ) \ 
v(q h )e^- e ^ 
u{q h )e-"t>* 
\u{q h )e 1 ^- 9h -' l >^ ) 



,—igh cos o h x 



(25) 



where (f>L,R is the superconducting phase on the left/right side 
of the normal region, associated with the broken U(l) sym- 
metry in the superconducting state. The macroscopic phase 
difference is defined as A<fi = <j)R — <pL- 

In order to proceed with an analytical treatment, we observe 
that g e /j = /is under the assumption that /is ^> {Ajj, At}. 
Moreover, the direction of the momentum that enters the argu- 
ment of the coherence functions {u, v} in the wavefunctions 
of Eqs. ( 23p5 i is of no significance since only the absolute 
value of the momentum enters At- The directional depen- 
dence has been separated out into the e 10 factors of the off- 
diagonal elements in Eq. ( 17 1. Thus, the problem effectively 
becomes equivalent to that of a conventional s-wave super- 
conductor with gap Ao = \Ajj — 3/isA t a/2|. So, while we 
have included both onsite and nearest neighbor interactions, 
thus giving rise to a combination of s-wave and p + lp-wave 
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pairing, the results would have been identical had we cho- 
sen only on-site or nearest neighbor interaction, as long as the 
nearest neighbor interaction gives rise to the extended s-wave 
symmetry in Eq. ff/J. This conclusion is supported by the find- 
ings of Jiang et al, who found that the Andreev conductance 
of a SN junction was identical for s-wave pairing and the ex- 
tended s-wave bond pairing (see their Fig. 3c). 

We now obtain the following energies for the Andreev 
bound states in the normal region: 



£± = ±\Au 
C(7) 



3/z s A t a/2|Vl-C(7)sin 2 (A0/2), 

P = ftsL/vp. (26) 



cos 2 7 



1 — sin 7 cos 2 (P cos ( 



Note that 9 = if A/i > fi s and 9 = 7 if A/i = 0, where 
7 = 9 e is the angle of incidence of quasiparticles. In the for- 
mer case, we regain the results of Maiti and Sengupta 1 ^ who 
studied the case of a thin and very strong barrier separating the 
two superconducting regions. In the case A/i = 0, we have 
no FVM between the graphene regions (/is = /xn = /i), and 
the results change accordingly. In this case the normalized 
Josephson current at zero temperature becomes 

r /2 d7C(7)cos 7 sinA< ? !) 
iJlh) = I 1 =) (2/) 



" /2 ^/l-C(7)sin 2 (A^/2) 

where Iq — bAq/K. In the following, we will study the 
dependence of the Josephson current on the (2) phase differ- 
ence, (if) doping level and length of the junction, and (2'2'z) 
the temperature-dependence. Since we consider the regime 
/i ^> {Ajj, At}, which means that fj, typically could lie in 
the range 10 — 100 meV, it is reasonable to expect that the 
role played by charge inhomogeneities such as electron-hole 
puddles may be disregarded in the main approximation. To be 
more specific, the local variations of the Fermi level 8/i should 
be of little importance as long as they satisfy 5fi <C /i. Exper- 
imentally, one has estimated^ 5/j, ~ 5 meV, which places a 
restriction on the appropriate values of /i. 

We now proceed to investigate the length-dependence of 
the critical current in Fig. |3ja). In all the plots presented in 
this section, the zero-temperature limit is assumed unless ex- 
plicitly stated otherwise. To operate within a valid regime of 
parameters, we restrict our attention to L/£ < 0.15, where 
£ = Hvf/Aq is the superconducting coherence length. As 
seen, the critical current displays an oscillating decay. This 
is a qualitatively new feature as compared to the results of 
Refs. [8 1 and |14|. For a strongly doped graphene junc- 
tion, but still with a FVM at the interfaces, Black-Schaffer et 
alP^ obtained numerically the length-dependence of the criti- 
cal current. From Fig. 5(b) in that work, one may see a hint of 
oscillations in the numerical data and it is clear that the cur- 
rent does not decay monotonously. The same data, but with 
more data points, are reproduced in Fig [TTJb) and from this 
and Fig. [4] it is clear that both the numerical and analytical 
data show an oscillatory dependence on L in the strong dop- 
ing regime including the case of no FVM. 

In Fig. [3jb), we give the current-phase relationship for 
several values of the length of the normal-region, using 



n/Ao = 50, in the case when there is no FVM. As seen, 
the current-phase relationship deviates slightly from the usual 
sinusoidal form with the phase difference providing the crit- 
ical current occurring at A<fi c 6 [tt/2,tt]. Qualitatively, 
Fig. [3jb) is in agreement with the numerical results shown 
in Fig. 4 of Ref. lfT2l . Note that by lowering /i and L, the 
current-phase relationship tends towards the functional form 
sin(A(/>/2)sgn{cos(A0/2)}. 

The analytical treatment above is valid for either A/i = 
(corresponding to zero FVM) or A/i ^> /is (corresponding to 
a barrier induced in the normal graphene region). Next, we 
consider a situation where the superconducting regions are 
strongly doped, while the normal region is only weakly or 
moderately doped. Since we have shown above that the effec- 
tive s+(p + lp) wave pairing near the Dirac points is formally 






-0.01 




-0.05 




-0.09 




-0.13 


§ \ 






1 






(b) \;l 




) 2 4 


6 



FIG. 3: (Color online) Josephson current for the s-wave symmetry 
with no FVM between the S and N regions (/is = /iN = M)- ( a ) 
Plot of the length-dependence of the critical current, (b) Plot of the 
current-phase relationship for /i/Ao = 50. 
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FIG. 4: (Color online) Josephson current for the s-wave symmetry 
with a strong FVM between the S and N regions (/ijv = /i). (a) 
Plot of the length-dependence of the critical current, (b) Plot of the 
current-phase relationship for L/£ = 0.1. 
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equivalent to the isotropic s-wave case, we can here use 8 : 



Cn sin A(f> 



h n=o ,/l-C n sin 2 (A0/2)' 



( n = kl[kl cos 2 (k n L) + sin 2 (k n L)} 1 , 

k n = q n = (n+ 1/2)tt/W. (28) 

Above, W denotes the width of the graphene junction, which 
we assume satisfies W ^S> L. We assume that the supercon- 
ducting regions are heavily doped, such that /is ^S> hn- In 
this case, the number of propagating modes in the supercon- 
ducting region is TV = /igW/w. We now proceed to investi- 
gate how the critical current depends on the length L of the 
junction. We choose W/£ = 30 and /j,g/Ao = 150. The re- 
sult is shown in Fig.Qa). In general, the actual magnitude of 
the current decreases with decreasing /i N since there are fewer 
propagating modes available when /in — * 0. At weak doping 
Hn I Ao ~ 1 — 10, the current decays as 1.33IqW/(ttL), in 
agreement with the findings of Ref. |8|. However, as the dop- 
ing level increases, oscillations appear as a function of L. We 
also give the current-phase relation for the case with a strong 
FVM in Fig. gb), using L/£ = 0.1. 



d x 2_ y 2- and d xy -wave 



As seen from Eqs. (FT), the situation becomes quite compli- 
cated when considering the d-wave symmetries that may arise 
upon considering nearest-neighbor pairing. At least, this is so 
when considering the DBdG equation in an atom-basis pic- 
ture, as has been done in the literature up to now. However, 
by transforming the DBdG-equation to the band-structure ba- 
sis (conduction and valence ir -bands), i.e. where the kinetic 
energy is diagonal, we are able to treat the d-wave symmetries 
analytically. This stems from the fact that the o?-wave symme- 
tries in this basis have a simply four-fold symmetry and their 
low energy expansions near the Dirac points result in even 
simpler p-waves. Below, we sketch this transformation. 

Our starting point is Eq. (7) in Ref. 11271 . where the super- 
conducting pairing is written in the band-picture as follows: 

H = Yl^ te i ~ V)c\ a c q(j + (-te q - fi)d q(T d qa ] 

qcr 

+ isin(q • aj - <p q ){c\^_ qi - 4 T cL gi )+H.c] (29) 
Here, we have defined 



= e n 



(30) 



Introducing a basis in the band-picture ^ 



[ C qp d q p C-qJ ■ ' 



d- q \], we may write the Hamiltonian as 



H = E q ^ g N q ^ q with 



fte q -n 



N„ 



-C* 



o -c q 

-te q - fj, -iS q 
iS* 

c* 



Co 



-te q + \l 
te a + 



iS q \ 



v 

where we have defined 

a 

S q = ^2 A j sin (<? • a j ~ ¥>«)) 



(31) 



(32) 



Note that N q is Hermitian, so we know that its eigenvalues 
will be real. As seen, the quantities {C qi S q } play the roles 
of the superconducting gaps for intra- and inter-band pairing, 
respectively. We are interested in the behaviour of N q near 
the Dirac points, effectively at the wavevector q = K± + k, 
where again |fe| < \K±\, and K ± = ±[0, 47r/(3v / 3a)]. By 
inserting this into N q and linearizing in k, we obtain for the 
d xy -wa\e symmetry featuring A = A t (0, 1, — 1) that 



C(K ± + k) = ±^^ 1 S(K ± +k) = 



A ^ h - (33) 



Effectively, this is a p x -wave pairing for Ck and p y -wave pair- 
ing for Sk- In an equivalent manner, we obtain for the d x 2_ y 2- 
wave case of A = A 4 (2, — 1, —1) that 

C(K ± + k) = T 3 -^, S(K ± + k) = -^-. (34) 

Let us also briefly mention in passing that for the extended 
s-wave case, A = A t (l, 1, 1), we obtain 



C{K ± + k) 



3A t \k\a 



S(K± + k) = 0, (35) 



which is consistent with the result of the previous section, 
i.e. the extended s-wave case gives rise to an effective conven- 
tional, fully gapped s-wave pairing near the Dirac points with 
a doping-dependent magnitude of the gap. The q-dependence 
in the entire BZ is shown for the c?-wave gaps in Fig. [5] 

The question is now: have we managed to simplify the ex- 
pressions compared to the atom-basis picture such that an ana- 
lytical approach has been rendered viable in the <i-wave case? 
At first sight, it appears that the situation is still rather compli- 
cated as there are two "gaps" in Eq. pi} , namely Ck and Sk- 
However, upon diagonalizing Nk to obtain its eigenvalues, we 
find 



E k = a\[tek + f5\KlJ 2 + S 2 )Y + C i 



(36) 



where again a, = ±1 refers to e-like and h-like particles, while 
f3 = ±1 refers to the conduction and valence band. From 



Eq. (36 1, it is clear that Sk simply renormalizes the chemical 
potential /i while Ck is the true superconducting gap. Assum- 
ing here a doped situation where /i > A t , we may certainly 
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d X y 




FIG. 5: (Color online) Plot of the (/-dependence of the two d- 
wave symmetry order parameters (C q ). Red color represent posi- 
tive sign, blue color negative sign, with zero being white. The first 
Brillouin zone is marked with a black line. Near the Dirac points, 
(<lx,qy) = (0, ±j^^), at the zone corners, effective p y - or p x - wave 
symmetries emerge in the low-energy regime. 



neglect Sk, compared to /x, and we therefore set Sk = in 
Eq. ( 3 1 1. The situation has now been considerably simplified. 



We are left with a single-gap superconductor with normal dis- 
persion <£fc and gap Ck where the gap has a simple p x - or 
Pj,-wave symmetry, which allows us to continue analytically. 

For the p^-wave symmetry, the situation becomes qualita- 
tively different from the conventional s-wave case since there 
are now zero-energy states (ZES) located near the interfaces. 
Since we are considering transport along the x-direction, the 
criterion for the existence of ZES is that the order parameter 
satisfies A(9) = — A(ir — 0). Clearly, this is the case for a 
p^-wave symmetry with A(9) ~ cos 9, whereas the p y -wavs 
case does not host any ZES. Let us here consider the case of 
no FVM. In the absence of ZES, the Andreev-bound states in 
the normal region may be written quite generally as 



±\A(i)\y[l -C(7)sin 2 (A0/2), 



(37) 



where again £(7) is the angularly resolved transmission co- 
efficient in the normal-state and Acf> is the phase-difference 
between the superconducting regions. It is seen that Eq. p7] i 
is formally equivalent to the Andreev bound-state spectrum 
for a conventional s-wave symmetry, with the only difference 
that the gap now has an angular dependence, i.e. A = A(7). 
As a result, one would expect qualitatively the same results for 
the Josephson current when comparing the p y -wave case with 
the s-wave case. Quantitatively, the magnitude of the current 
would be reduced due to angular averaging over the gap. 

We now consider the p^-wave symmetry featuring ZES, 
and one then in general findP 



then becomes 

Ij/I = 2 l 7sin(A0/2)sgn{cos(A</>/2)}, 

J= d 7 cos 2 7VC(7)) 

J-tt/2 



C(7) = 



cos 2 7 



1 — sin 7 cos 2 (P cos 7) 



, P = [j,L/vf- (39) 



We proceed to discuss the length-dependence and phase- 
dependence of this Josepshon current, and especially inves- 
tigate how it differs from the conventional s-wave case where 
there are no ZES. As seen in Fig. |6|a), the current still dis- 
plays oscillations as a function of the width L in the strongly 
doped case 11 ^S> Ao, but the oscillation-amplitude is consid- 
erably smaller than for the s-wave symmetries. In Fig. |6jb), 
we give the current-phase relationship for the Josephson junc- 
tion. As seen, there is now an abrupt crossover at A<fi = it 
which should be contrasted with the smooth behavior in the 
s-wave case shown in Fig. |3jb). An interesting aspect is that 
in the present case of d-wave symmetry, the junction energy 
is minimized at A<fi = it, while in the s-wave case the mini- 
mum of free energy occurs at a phase difference </>o which lies 
between and ir. 
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FIG. 6: (Color online) Josephson current for the d^-wave symme- 
try (giving rise to an effective p^-wave symmetry at the Dirac points) 
with no FVM between the S and N regions, (a) Plot of the length- 
dependence of the critical current, (b) Plot of the current-phase rela- 
tionship for h/Aq — 50. 

We proceed to investigate the Josephson current when there 
is a substantial FVM between the S and N regions. In this 
case, we obtain from Eq. (|38]> that 



A N 

71 = 



C n sin(A(/./2)sgn{cos(A< ? !>/2)}, 



(40) 



±|A( 7 )|VC(7)cos(A0/2), 



(38) 



where, A (7) = Aq cos 7. The normalized Josephson current 



where t n models the angular dependence of the gap. For the 
most interesting p^-wave case, we choose t n — cos(^). The 
resulting critical current is shown in Fig. |7ja), using the pa- 
rameters hs/Aq — 150 and W/£ = 30. As seen, the current 
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FIG. 7: (Color online) Josephson current for the d xy -wave symmetry 
(giving rise to an effective p^-wave symmetry at the Dirac points) 
with a strong FVM between the S and N regions, (a) Plot of the 
length-dependence of the critical current, (b) Plot of the current- 
phase relationship for L/£ = 0.1. 



decays like L^ 1 in the regime fi^ / Aq ~ 1 — 10 while os- 
cillations appear upon increasing /ijy further, just as in the 
s-wave case. In agreement with Fig. |6ja), it is seen that the 
oscillation-amplitude is smaller than in the s-wave case. The 
current-phase relationship in the case of a strong FVM for the 
e?-wave symmetries is shown in Fig. [7jb). 
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B. Numerical results 

In this section we will report on the self-consistent numeri- 
cal results for SNS graphene junctions with both conventional 
s-wave and unconventional SB pairing contacts. Our main 
focus is to complement the analytical work and point to sit- 
uations where a self-consistent approach is a must in order 
to capture the correct behavior. We are also able to extract 
additional data not available through an analytical approach, 
such as the proximity effect and LDOS, and we will start with 
reporting these below. This will provide the necessary back- 
ground to interpret the Josephson current and how it in some 
cases differs significantly from the analytical result. 

As shown above through analytical work, the extended s- 
wave solution for SB pairing contacts manifests itself as an 
on-site pairing gap at non-zero doping, which is always the 
case in the S regions. We therefore choose to not focus on this 
solution and only study the two distinct d-wave solutions in 
Eq. |2]i and compare these with on-site s-wave contacts. More 
detailed results for the on-site s-wave solution can be found 
in earlier work by some of the authors.^ 

Figure|8]shows the proximity effect in terms of the normal- 
ized pair amplitude for several superconducting symmetries 
and interface combinations and at multiple doping levels in 
the N region. Solid black lines are the results for on-site s- 
wave pairing which show a pronounced depletion of pairs on 
the S side of the interfaces and the accompanied leakage of 
pairs into the N regions. This depletion/leakage is stronger 



FIG. 8: (Color online) Proximity effect in terms of normalized pair 
amplitudes when ^i(N) = eV (a), /Lt(N) = 0.7 eV (b), and jtt(N) = 
1.5 eV (no FVM) (c). Conventional s-wave contacts (solid black), 
d a .2_ J( 2-wave contacts with zigzag interfaces (dashed red), d^-wave 
contacts with zigzag interfaces (dotted blue), and d a! 2_ !/ 2-wave con- 
tacts with armchair interfaces (dash-dotted green). The width of the 
N region is L — 0.42£ and the interfaces are marked with vertical 
black lines. 



at higher doping levels. The red dashed curve is the results 
for the zigzag interface with d x 2_ y 2-wave symmetry contacts, 
i.e. an effective p y -wave symmetry at low energies, and it 
displays a much weaker proximity effect compared to the s- 
wave contacts. In fact, at zero doping in N there is only a 
very small, and interestingly, completely flat, non-distance 
dependent, superconducting pair amplitude inside the junc- 
tion. The small oscillatory pattern present, especially at lower 
doping, for this junction can be attributed to charge fluctua- 
tions due to the FVM at the interface and is seen also for the 
s-wave solution, 12 though it is less pronounced there. The 
green dash-dotted and blue dotted lines are the results for the 
d x 2_ y 2-wa\e pairing on the armchair interface and the dxy- 
wave contacts on the zigzag interface, respectively. Both of 
these have an effective p^-wave symmetry at low energies and 
should therefore to a good approximation be the same, as is 
also seen here. We will thus from hereon only report results 
for the armchair (4.2 _„2 -wave pairing. The depletion of pair 
amplitude on the S side of the interface is notably larger for 
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these contacts compared to s-wave contacts, but the induced 
pair amplitude in N is nevertheless not enlarged. Colloquially 
speaking this means that these contacts have in total lost more 
pairs than junctions with conventional contacts. This effect 
can be attributed to the formation of quasiparticle ZES at the 
SN interfaces. 

The presence of ZES for the d^-wave symmetry on the 
zigzag interface, or equivalently the d x 2_ J/ 2-wave symmetry 
on the armchair interface, is the most prominent differences 
between the two different <i-wave symmetries in the analyti- 
cal framework. The TB BdG framework allows for a direct 
access to the LDOS and therefore any ZES formation. Fig- 
ure [9] shows the results for s-wave contacts (a, d) and d x i_ y i- 
wave contacts on the zigzag (b, e) and armchair interfaces (c, 
f). For conventional s-wave contacts we see, as expected, a 
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FIG. 9: LDOS plots for conventional s-wave contacts (a, d), d x 2_ y 2- 
wave contacts with zigzag interfaces (b, e) and armchair interfaces (c, 
f). The doping level in the N region is zero (a, b, c) and moderately 
doped at 0.7 eV (d, e, f). The energy scale has been normalized by 
the order parameter Alt. Black color represent 2 states/eV/unit cell. 

full gap in the S regions and for short enough junctions this 
gap persists inside N for both zero (a) and moderately doped 
(d) N regions. However, the situation is quite different for the 
d-wave contacts. Here, the order parameter in the band struc- 



ture basis, C q in Eq. (32 1, has nodes on the Fermi surface 
leading to a familiar V-shaped DOS in the gap. In addition 



to this feature, we also see pronounced peaks in the LDOS 
at zero energy at the SN armchair interfaces for £^2^2 -wave 
contacts as predicted. These ZES are most prominent in the 
large FVM limit (c) but exists even when the N region is mod- 
erately doped and the FVM is smaller (f). At no FVM there is 
no trace of ZES at the interface. This case corresponds to a di- 
minishing interface barrier Z and even in regular BTK-theory 
for d-wave SN junctions there is then no distinct signatures of 
the ZES since the transmission probability is unity for Z = 
(see e.g. Ref. Il34l ). Note that the lighter color throughout the 
N region in all junctions simply reflect the fact that the DOS 
is lower here since /ijv < [is- 

In Figs. [8] and [9] we have seen that different order param- 
eter symmetries in the contacts alter the proximity effect and 
LDOS significantly and it is therefore expected that there will 
be a large difference in the Josephson currents as well. Here 
we will focus on two cases; the dependence of the Joseph- 
son current on the phase difference across the junction and 
on junction length. Figure 10 shows the dependence of the 
Josephson current / on the phase difference Acf> across the 
junction for different doping levels in N for a junction of 
length L = 0.42£. Here Atfi is the phase difference just across 
the N region which will always be smaller than the largest 
phase difference between the two contacts, tt, for any finite 
Josephson current. This is why the numerical results, espe- 
cially for larger currents, cannot be extended to large Acf>- 
values. The numerical results are explicitly calculated for in- 
finitely wide junctions, but by letting W — > 00 for the analyt- 
ical results we are able to compare analytical and numerical 
solutions in the zero and moderately doped cases. In the case 
of no FVM, the analytical result does not depend on the width 
and a comparison is much less straightforward and therefore 
not plotted in Fig.[T0|c). 

Let us first focus on the numerical results alone, which are 
the solid lines in Fig. 10 For zero doping in N, d x 2_ y 2-wscve 
contacts on the armchair interface (green,+) have the highest 
current. This is expected since this is the only configuration 
with ZES at the SN interfaces and these states will strongly 
intensify the tunneling current through the junction. 35 At least 
in terms of the Josephson current, the formation of the ZES 
clearly makes up for the relative loss of pair amplitude seen 
in Fig. [8] We note that the relative enhancement in the current 
for d 2 ,2_ y 2-wave contacts on armchair interface compared to 
s-wave contacts is reduced when the FVM at the interface is 
reduced. This is also to be expected as the strength of the ZES 
diminishes with increased doping level in N. In fact, when the 
ZES disappear at no FVM, the Josephson current is lower than 
for the similarly strong s-wave pairing contacts. Of all the 
symmetries investigated, dj.2_.y2 -wave pairing on the zigzag 
interface (red, o) has the lowest current in all three doping 
regimes. This is consistent with its smallest proximity effect 
as seen in Fig. [8] Finally, also note that the current for all three 
different pairings increases with increasing doping level in N. 

Let us now comment on the explicit Atfi dependence on 
the current and also make a comparison with the analytical 
results. The s-wave contacts show a distinct non-sinusoidal 
dependence with the numerical results closely tracking the 
analytical results at low doping levels. At moderately dop- 
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FIG. 10: (Color online) I vs. A(j> for ^ N = eV (a), 0.7 eV (b), and 
1.5 eV (no FVM) (c) for conventional s-wave contacts (black, x), 
d a .2_j / 2-wave contacts with zigzag interfaces (red, o) and armchair 
interfaces (green, +). The current is given in units of 7o = eW/ (h£) 
when vf is set to 1. Analytical results in (a, b) are shown with dashed 
lines and calculated for the width W = 30£, which approximates the 
infinite width of the numerical results. 



ing the analytical results peak at a higher A0-value than the 
self-consistent numerical results, around 0.727T compared to 
0.567T. However, the discrepancy is only moderate, both in 
terms of critical A<fi and I, as already reported in Ref. Ifl2l . 
For the <i-wave contacts the situation is, however, quite differ- 
ent. Especially big is the discrepancy for the case with ZES 
at the interfaces. Here, the analytical results show a distinct 
sgn{cos(A0/2)} sin(A0/2) behavior, which is a direct con- 
sequence of the presence of the ZES, and a very large criti- 
cal current for all doping levels. The self-consistent results 
instead display an almost sin(A^)-like curve with a signifi- 
cantly lower critical current. We believe this large difference 



can be attributed to the importance and strength of the ZES. 
While the ZES are present even in the self-consistent solu- 
tion and there causes an enhancement of the Josephson cur- 
rent, its importance seem to be significantly reduced within a 
self-consistent approach. This would explain both the smaller 
amplitude and the more traditional sin(A</>)-curve. Since the 
effect of the ZES is decreased with increasing doping level 
the discrepancy between the analytical solution and the self- 
consistent work slightly decreases in terms of maximum cur- 
rent, but the different A0-dependence will remain as long as 
the ZES are non-zero in either approach. It might be worth 
mentioning here that ZES have been found to be very sensitive 
to interface properties. Similar numerical solutions schemes 
as applied here for <i-wave contacts on the square lattice have 
found that the ZES can quickly diminish in the presence of 
random interface potentials, simulating disorder,^ or for non- 
straight interfaces destructive interference between different 
lattice sites can completely kill the ZESP. It is therefore very 
likely that the reduced importance of the ZES in our numerical 
solution stems from the fact that even a completely smooth, 
disorder free armchair interface as present in our calculations, 
is not absolutely flat in the continuum-sense. 

However, also the d :E 2„ y 2-wave symmetry on the zigzag in- 
terface shows pronounced differences between the analytical 
and self-consistent solutions. Here, the phase dependence is 
similar, with the critical phase difference slightly increasing 
with doping level in N, but the current for these junctions is 
significantly enhanced in the self-consistent solution. At the 
Dirac point there is a ~ 60% increase in the critical current 
compared to the analytical solution. This enhancement in- 
creases with doping level in N, being ~ 80% at 0.7 eV doping. 
It is interesting that we see this enhancement in the current for 
the self-consistent solution despite the proximity effect caus- 
ing less leakage into N for these contacts compared to s-wave 
pairing. That also means less depletion of the order param- 
eter in the S region of the interface which will give rise to a 
stronger tunneling and, apparently, this effectively even over- 
compensates for the relative lack of pair amplitude in N. 

The above shows that for unconventional contacts a self- 
consistent approach is necessary in order to accurately deter- 
mine the Josephson current, both in terms of phase depen- 
dence and critical current. We have identified at least two 
causes for this. First, at SN interfaces where ZES are present, 
a self-consistent calculation is necessary in order to properly 
evaluate the importance of the ZES to the overall transmis- 
sion of the junction. Second, the proximity effect itself can 
be considerably different for unconventional contacts as com- 
pared to s-wave contacts. This gives rise to the main source 
of discrepancy in the results between the analytical and self- 
consistent approach in the case of the d x 2_ y 2-wa\e symmetry 
on the zigzag interface. 

Next, we report on the junction length dependence of the 
Josephson current. The analytical approach used here is only 
reliable in the short junction regime where the current is 
mainly carried by the Andreev bound states. For the self- 
consistent TB BdG approach such a limitation does not exist 
and we can therefore extend the results to much longer junc- 
tions. On the other hand, limited computational time requires 
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relatively short superconducting coherence lengths and with 
the parameters used here the shortest junctions we can study 
have L — 0.04£. Also, for these short junctions there will be 
a substantial current across the junction and, as a result, the 
available phase difference A(j> across the junction is drasti- 
cally reduced making it hard to reach A<fr c for very short junc- 
tions. This means that we effectively do not have any results 
for junctions shorter than ~ 0.1£ for [i — eV and ~ 0.25£ 
for [i = 0.7 eV. For relatively large junction lengths the ana- 
lytical results unfortunately start to become questionable and 
we can therefore not make a comparison between the analyt- 
ical and the self-consistent results over any extended length 
scale. Figure 11 shows the critical current I c vs. length of 
junction L for all three cases of symmetries of the pairing at 
doping levels eV (a) and 0.7 eV (b) in N. At zero doping in N 
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FIG. 11: (Color online) I c vs. L for /i N = eV (a) and 0.7 eV (b) 
for conventional s-wave contacts (black, x), d a .2_ H 2-wave contacts 
with zigzag interfaces (red, o) and armchair interfaces (green, +). 
The current is given in units of Jo = eW/ (h£) when vf is set to 
1 . Analytical results are shown with dashed lines and calculated for 
the width W = 30£, which approximates the infinite width of the 
numerical results. 



(a) s-wave contacts show a relatively good agreement between 
the analytical and self-consistent results, although the length 
dependence is somewhat stronger in the self-consistent case. 12 
A simple least square fit to the functional form I c = CL^ 13 
gives (3 — 1.3 for the self-consistent curve but only 0.8 for the 
analytical results. The d 2 .2_ y 2-wave contacts on the zigzag in- 
terface show a similar increased length dependence for short 
junctions with (3 = 2 for the self-consistent results and 1.3 
for the analytical counterpart. The d x 2_ y 2-wave contacts on 
the armchair interface have the weakest L-dependence with 
only j3 = 1.1 for the self-consistent results. The analytical 
results in this case give a bad fit to a power-law dependence as 
the critical current levels off for short junctions. There are no 
noticeable oscillations in the critical current as a function of 
junction length in the undoped regime for either of s-wave or 
the £^2^2 -wave contacts with the zigzag interfaces for either 
solution method. However, d x 2_ v 2-wave contacts with the 



armchair interface exhibit peaks in the critical current at every 
6 unit cells or 0.14£ in the self-consistent approach. These 
are most pronounced for short junctions, but also exist for 
the longest junctions we have studied. Since the analytical 
solution does not capture these oscillations, we attribute the 
peaks to resonance transmission due to changes in the relative 
strength or interactions of the ZES. 

At moderate doping and thus a smaller FVM (b) we see 
a fairly large discrepancy in the critical value of the current 
between the analytical and the self-consistent solutions for 
all three contacts, as also seen in Fig. 10 The length de- 



pendence is however very similar in the two solutions for 
both s-wave and d x 2_ J( 2-wave contacts on the zigzag inte- 
face. Here (3 = 0.5 and 0.85, respectively. For the junc- 
tion with ZES, the length dependence is harder to estimate 
because of the lack of short junction self-consistent data but 
the results approximately vary from (3 — 0.2 to 0.5 when 
self-consistency is included. At moderate doping both the an- 
alytical and the self-consistent data show long wave length 
oscillations for especially the s- and c? 2 .2_ !/ 2-wave contacts 
with the zigzag interfaces. While it is not straightforward 
to compare frequencies and amplitudes of the analytical and 
self-consistent approaches there seem to be an overall agree- 



ment. For cL 



2 -wave contacts on the armchair interface, we 



again find short wave lengths resonance transmission peaks in 
the self-consistent solution, but we can in this case not easily 
identify any distinct frequency. 

Not only do the I vs. A<fi and I c vs. L curves differ most 
for the 6^2 _ y 2 -wave contacts on the armchair interface, but 
Fig. 12 shows that there is a strong junction length dependence 



on the critical phase difference A<j) c where the maximum cur- 
rent is reached. As seen in Figs. 6]|7 A<p c = tt for all junction 




FIG. 12: (Color online) I c (black, left axis) and A(p c (red, right axis) 
vs. L for /xn = eV (a) and 0.7 eV (b) for the d x 2_ !/ 2-wave solution 
on the armchair interface. 

lengths in the analytical solution in this case, whereas the self- 
consistent solution shows A<fi c increasing from around 0.37T 
in the shortest junctions we could study to above 0.67T in the 
longest junctions. The spread in A<fr c could be even larger, 
but the curves seem to level out at large junction lengths. In- 
terestingly, this does not only mean that the maximum current 
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is not reach at a 7r phase difference, but also that the phase 
difference is advanced compared to the conventional sin(A0) 
curve for short junctions. It might be worth noting that there is 



a fair amount of fluctuations in the data in Fig. 12 Part of this 



is due to the way we measure the Josephson current. We apply 
a fixed phase difference across the whole SNS structure, thus 
technically we are only able to give an interval for the critical 
current and phase difference. However, the peaks in critical 
current seen in undoped junctions for every 6 unit cell are not 
a numerical artifact, and as seen, they are also not correlated 
with the wiggles in the critical phase difference. We admit 
that it is somewhat harder to draw the same conclusion in the 
moderately doped case. There is no similar junction length 
dependence for s-wave or d x 2_j,2-wave pairing on the zigzag 
interface, but here A0 C is fixed over the whole range of L we 
investigate. 



IV. CONCLUSIONS 

In summary, we have studied the Josephson current in a 
graphene superconductor/normal/superconductor (SNS) junc- 
tion. Superconductivity is induced in graphene by means 
of the proximity effect from a superconducting host mate- 
rial. In particular, we assume that superconducting con- 
tacts are deposited on top of the graphene sheet as re- 
alized experimentally. 4 Whereas previous works on the 
Josephson current in graphene have mainly treated only an 
isotropic s-wave symmetry for the superconducting order 
parameter, ^ 8 ! 9 ! 10 ! 11 ! 12 ! 13 ! 14 ! 16 ! 18 ! 20 ! we here also present an anal- 
ysis for anisotropic pairing that arises due to nearest-neighbor 
interactions. This latter pairing can either have an extended 
s-wave symmetry or belong to any linear combination of d xy 
and d x 2_ y 2. 

We show that a junction with extended s-wave symmetry, 
which displays an effective p x + ip y -wave symmetry near the 
Dirac points, is equivalent to a junction with on-site, isotropic 
s-wave pairing. While s-wave pairing has been studied be- 
fore, we here report on newly found oscillations in the critical 
current as function of junction length in both junctions with 
no Fermi vector mis-match (FVM) at the SN interfaces and in 
heavily doped junctions with a strong FVM. 

For the case of d-wave superconducting contacts we limit 
our investigation to considering only d x 2_ y 2-wa\e pairing on 
the zigzag and armchair interfaces. Since all pairings are 
induced from on-top deposited contacts, different symmetry 
choices simply correspond to different orientations of the con- 
tacts relative to the graphene sheet. These chosen symmetries 
give rise to effective p x - and p y -wa\e pairing, respectively, at 
low energies. Therefore d x 2_ y 2-wa\e pairing on the armchair 



interface is at low energies equivalent to d^-wave pairing on 
the zigzag interface and so on, making our study quite general. 
In an experimental setup any chiral interface on the graphene 
could be realized. However, we believe that this situation can 
at least qualitatively be predicted from our results and will 
mainly depend on the presence of zero energy states (ZES) 
at the interfaces. These states appear if the order parameter 
changes sign when the angle of incidence for the quasiparti- 
cles on the SN interface is changed from 8 to tt — 8 as is the 
case for d x 2_ y 2-wave pairing on the armchair interface. 

We calculate the Josephson current both analytically and 
numerically by a self-consistent approach for all the above 
symmetries. Whereas there is good agreement between the 
two treatments for the s-wave superconducting order parame- 
ters, there is a pronounced deviation between the two methods 
for anisotropic pairing, in particular, when ZES are present. 
These states at zero energy will easily dominate the trans- 
port through the junction, and while they are present even in 
the numerical results, their effect on the Josephson current is 
strongly suppressed when self-consistency is achieved. One 
easily identified source for this deviation is the first order ex- 
pansion to p-wave symmetries around the Dirac points done 
in the analytical treatment. As seen in Fig. [5] the d-wave or- 
der parameters differ from pure p- waves at higher energies. 
Since the contacts are likely to induce a rather heavy doping 
into the graphene in the S regions one might argue that this 
effect is not negligible. However, as seen in Fig. [8] the prox- 
imity effect is remarkably similar for the d^-wave on zigzag 
interface and the d x 2_ y 2-wa\e on the armchair interface, thus 
pointing to the fact that this is not a major source of deviation 
between the analytical and the self-consistent treatment. It is 
instead the self-consistency for the proximity effect inside the 
junction that is the crucial component. Therefore, a numeri- 
cal, self-consistent calculation is required in order to properly 
address the transport properties of graphene when the super- 
conducting pairing is anisotropic in fc-space, especially when 
zero-energy states are present at the interface. 
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